##################################################
## Author: Stephanie L. DeMora & Benjamin Newman
## Project: Adolescent Exposure to Economic Inequality and Belief in the "American dream" in Entering Adulthood
## Purpose: Replication Code
## Date: 2025-11-20
##################################################

pkgs <- c("tidyverse", "haven", "janitor", "interflex", "lmerTest", "DT", "here", "jtools",
          "ggeffects", "ggpubr", "clubSandwich", "stargazer", "sandwich", "ggplot2", "zipcodeR", "effectsize")
lapply(pkgs, library, character.only = TRUE)

df1 <- readRDS("Data_Pt1.rds")
df2 <- readRDS("Data_Pt2.rds")
df <- rbind(df1, df2); rm(df1, df2); gc()

# ---------------------------------------------------------------------------------------------------------------
# ---------------------------------------- APPENDIX A -----------------------------------------------------------
# ---------------------------------------------------------------------------------------------------------------
# See "Replication Code - SI Appendix A Only.R, as this section uses additional data sources.

# ---------------------------------------------------------------------------------------------------------------
# ---------------------------------------- APPENDIX B -----------------------------------------------------------
# ---------------------------------------------------------------------------------------------------------------

#-----------------------------------------------------------------------------
# Table S4. Descriptive Statistics of TFS between 2005 and 2008
#-----------------------------------------------------------------------------
df %>%
  filter(YEAR == 2005 | YEAR == 2006 | YEAR == 2007 | YEAR == 2008) %>%
  get_summary_stats(Belief_in_Meritocracy, Gini, 
                    `Zip Education`, `Zip Income`, `Zip Unemployment`, `Zip White`, `Zip Density`,
                    Female, White, Citizen, Political_Orientation, `No Religious Affiliation`, Attend_Church,
                    Income, `Mothers Education`, `Fathers Education`, `Conservative Parent Career`, `Percent Voted Bush`, 
                    show = c("mean", "sd"))

# ---------------------------------------------------------------------------------------------------------------
# ---------------------------------------- APPENDIX C -----------------------------------------------------------
# ---------------------------------------------------------------------------------------------------------------

#-----------------------------------------------------------------------------
# Table S5. Underlying models for the left (Model 1) and middle (Model 2) panel of Figure 1 
#-----------------------------------------------------------------------------
mod1 <- lmer(Belief_in_Meritocracy ~ Gini + `Zip Education` + `Zip Income` + `Zip Unemployment` + `Zip White` + `Zip Density` +
               Female +   White + Citizen + Political_Orientation + `No Religious Affiliation` + Attend_Church +
               Income + `Mothers Education` + `Fathers Education` + `Conservative Parent Career` + `Percent Voted Bush` + 
               (1 | YEAR) + (1 | School_Code), data = df)
mod2 <- lmer(Belief_in_Meritocracy ~ Ratio_8020 + `Zip Education` + `Zip Income` + `Zip Unemployment` + `Zip White` + `Zip Density` +
               Female +   White + Citizen + Political_Orientation + `No Religious Affiliation` + Attend_Church +
               Income + `Mothers Education` + `Fathers Education` + `Conservative Parent Career` + `Percent Voted Bush` + 
               (1 | YEAR) + (1 | School_Code), data = df)

export_summs(mod1, mod2, digits = 4, error_format = "({std.error})")

#-----------------------------------------------------------------------------
# Table S6. Underlying model for the right panel of Figure 1
#-----------------------------------------------------------------------------
df %>%
  mutate(Zip_Education = `Zip Education`,
    Zip_Income = `Zip Income`,
    Zip_Unemployment = `Zip Unemployment`,
    Zip_White = `Zip White`,
    Zip_Density = `Zip Density`,
    Percent_Voted_Bush = `Percent Voted Bush`,
    School_Code = `School Code`,
    No_Religion = `No Religious Affiliation`,
    Mothers_Education = `Mothers Education`,
    Fathers_Education = `Fathers Education`,
    Conservative_Parent = `Conservative Parent Career`) %>%
  select(YEAR, Income, School_Code, Belief_in_Meritocracy, Gini, 
         Zip_Education, Zip_Income, Zip_Unemployment, Zip_White, Zip_Density,
         Percent_Voted_Bush, No_Religion, 
         Mothers_Education, Fathers_Education, Conservative_Parent,
         Female, Age, White, Citizen, Political_Orientation, Attend_Church, pmore75zip, pless25zip) %>%
  filter(complete.cases(Belief_in_Meritocracy, Gini, pless25zip, pmore75zip)) %>%
  as.data.frame() -> interaction_data

Plot1C <- interflex(estimator = "binning",Y = "Belief_in_Meritocracy", D = "pless25zip", X = "pmore75zip", 
                   Z = c("Female", "White", "Citizen", "Income",
                         "Political_Orientation", "Attend_Church", "No_Religion",
                         "Mothers_Education", "Fathers_Education", "Conservative_Parent",
                         "Zip_Education", "Zip_Income", "Zip_Unemployment", "Zip_White",
                         "Zip_Density", "Percent_Voted_Bush"), cl = c("YEAR", "School_Code"),
                   data = interaction_data, na.rm = TRUE,
                   Ylabel = "", Dlabel = "", Xlabel="% above $75k",show.grid = F,
                   cex.main = 1.2,
                   theme.bw = TRUE, bin.labs = F,ylab = "Marginal effect of % below $25k",
                   parallel=TRUE)


round((summary(Plot1C$model.binning)$coefficients), digits = 4)

#-----------------------------------------------------------------------------
# Figure 1, Panel C
#-----------------------------------------------------------------------------
Plot1C

#-----------------------------------------------------------------------------
# Table S7. Underlying model for the left panel of Figure 2
#-----------------------------------------------------------------------------
df %>%
  mutate(Gini = Gini,
    Zip_Education = `Zip Education`,
    Zip_Income = `Zip Income`,
    Zip_Unemployment = `Zip Unemployment`,
    Zip_White = `Zip White`,
    Zip_Density = `Zip Density`,
    Percent_Voted_Bush = `Percent Voted Bush`,
    School_Code = `School Code`,
    No_Religion = `No Religious Affiliation`,
    Mothers_Education = `Mothers Education`,
    Fathers_Education = `Fathers Education`,
    Conservative_Parent = `Conservative Parent Career`) %>%
  select(YEAR, Income, School_Code, Belief_in_Meritocracy, Gini, 
         Zip_Education, Zip_Income, Zip_Unemployment, Zip_White, Zip_Density,
         Percent_Voted_Bush, No_Religion, 
         Mothers_Education, Fathers_Education, Conservative_Parent,
         Female, Age, White, Citizen, Political_Orientation, Attend_Church) %>%
  filter(complete.cases(Belief_in_Meritocracy, Gini, Income)) %>%
  as.data.frame() -> interaction_data

Plot2A <- interflex(estimator = "binning",Y = "Belief_in_Meritocracy", D = "Gini", X = "Income", 
                   Z = c("Female", "White", "Citizen",
                         "Political_Orientation", "Attend_Church", "No_Religion",
                         "Mothers_Education", "Fathers_Education", "Conservative_Parent",
                         "Zip_Education", "Zip_Income", "Zip_Unemployment", "Zip_White",
                         "Zip_Density", "Percent_Voted_Bush"), cl = c("YEAR", "School_Code"),
                   data = interaction_data, na.rm = TRUE, ylim = c(-0.4, 0.2),
                   Ylabel = "Belief in Meritocracy", Dlabel = "GINI", xlab = "",
                   cex.main = 1.2,
                   theme.bw = TRUE, show.grid = FALSE,
                   parallel=TRUE)

round((summary(Plot2A$model.binning)$coefficients), digits = 4)

#-----------------------------------------------------------------------------
# Figure 2, Panel A
#-----------------------------------------------------------------------------
Plot2A

#-----------------------------------------------------------------------------
# Table S8. Underlying model for the right panel of Figure 2
#-----------------------------------------------------------------------------
df %>%
  mutate(Gini = Gini,
    Zip_Education = `Zip Education`,
    Zip_Income = `Zip Income`,
    Zip_Unemployment = `Zip Unemployment`,
    Zip_White = `Zip White`,
    Zip_Density = `Zip Density`,
    Percent_Voted_Bush = `Percent Voted Bush`,
    School_Code = `School Code`,
    No_Religion = `No Religious Affiliation`,
    Mothers_Education = `Mothers Education`,
    Fathers_Education = `Fathers Education`,
    Conservative_Parent = `Conservative Parent Career`) %>%
  select(YEAR, Income, School_Code, Belief_in_Meritocracy, Ratio_8020, 
         Zip_Education, Zip_Income, Zip_Unemployment, Zip_White, Zip_Density,
         Percent_Voted_Bush, No_Religion, 
         Mothers_Education, Fathers_Education, Conservative_Parent,
         Female, Age, White, Citizen, Political_Orientation, Attend_Church, HOMEZIP) %>%
  filter(complete.cases(Belief_in_Meritocracy, Ratio_8020, Income)) %>%
  as.data.frame() -> interaction_data


Plot2B <- interflex(estimator = "binning",Y = "Belief_in_Meritocracy", D = "Ratio_8020", X = "Income", 
                   Z = c("Female", "White", "Citizen",
                         "Political_Orientation", "Attend_Church", "No_Religion",
                         "Mothers_Education", "Fathers_Education", "Conservative_Parent",
                         "Zip_Education", "Zip_Income", "Zip_Unemployment", "Zip_White",
                         "Zip_Density", "Percent_Voted_Bush"), cl = c("YEAR", "School_Code", "HOMEZIP"),
                   data = interaction_data, na.rm = TRUE,
                   Ylabel = "Belief in Meritocracy", Dlabel = "Ratio_8020", xlab = "",
                   cex.main = 1.2, 
                   theme.bw = TRUE, show.grid = FALSE,
                   parallel=TRUE)

round((summary(Plot2B$model.binning)$coefficients), digits = 4)

#-----------------------------------------------------------------------------
# Figure 2, Panel B
#-----------------------------------------------------------------------------
Plot2B

# ---------------------------------------------------------------------------------------------------------------
# ------------------------Appendix-D: Alternative main model specifications -------------------------------------
# ---------------------------------------------------------------------------------------------------------------

#-----------------------------------------------------------------------------
# Table S9.
#-----------------------------------------------------------------------------
mod1 <- lmer(Belief_in_Meritocracy ~ Gini + `Zip Education` + `Zip Income` + `Zip Unemployment` + `Zip White` + `Zip Density` +
               Female +   White + Citizen + Political_Orientation + `No Religious Affiliation` + Attend_Church +
               Income + `Mothers Education` + `Fathers Education` + `Conservative Parent Career` + `Percent Voted Bush` + 
               (1 | YEAR) + (1 | School_Code) + (1 | HOMEZIP), data = df)

mod2 <- lm(Belief_in_Meritocracy ~ Gini + `Zip Education` + `Zip Income` + `Zip Unemployment` + `Zip White` + `Zip Density` +
                      Female +   White + Citizen + Political_Orientation + `No Religious Affiliation` + Attend_Church +
                      Income + `Mothers Education` + `Fathers Education` + `Conservative Parent Career` + `Percent Voted Bush` + as_factor(YEAR) +
                      as_factor(School_Code), data = df)

mod3 <- lmer(Belief_in_Meritocracy ~ Gini + `Zip Education` + `Zip Income` + `Zip Unemployment` + `Zip White` + `Zip Density` +
                    Female +   White + Citizen + Political_Orientation + `No Religious Affiliation` + Attend_Church +
                    Income + `Mothers Education` + `Fathers Education` + `Conservative Parent Career` + `Percent Voted Bush` + 
                    (1 | YEAR) + (1 | School_Code) + (1 | COUNTY), data = df)

export_summs(mod1, mod2, mod3, digits = 4, error_format = "({std.error})")

#-----------------------------------------------------------------------------
# Table S10. 
#-----------------------------------------------------------------------------
mod1 <- lme4::lmer(Belief_in_Meritocracy ~ Gini + `Zip Education` + `Zip Income` + `Zip Unemployment` + `Zip White` + `Zip Density` +
               Female +   White + Citizen + Political_Orientation + `No Religious Affiliation` + Attend_Church +
               Income + `Mothers Education` + `Fathers Education` + `Conservative Parent Career` + `Percent Voted Bush` + 
               (1 | YEAR) + (1 | School_Code) + (1 | HOMEZIP), data = df)
se_mod1   <- sqrt(diag(vcov(mod1)))

mod2 <- lm(Belief_in_Meritocracy ~ Gini + `Zip Education` + `Zip Income` + `Zip Unemployment` + `Zip White` + `Zip Density` +
              Female +   White + Citizen + Political_Orientation + `No Religious Affiliation` + Attend_Church +
              Income + `Mothers Education` + `Fathers Education` + `Conservative Parent Career` + `Percent Voted Bush` + YEAR + as_factor(School_Code), data = df)
cl_se_mod2 <- sqrt(diag(vcovCL(mod2, cluster = ~HOMEZIP)))

stargazer(mod1, mod2,
          se = list(se_mod1, cl_se_mod2),
          no.space = TRUE,
          omit.stat = c("LL","ser","f"),
          column.sep.width = "50pt",
          align = TRUE,
          type = "text")

#-----------------------------------------------------------------------------
# Table S11. 
#-----------------------------------------------------------------------------
mod1 <- lm(Belief_in_Meritocracy ~ Gini, data = df)
mod2 <- lmer(Belief_in_Meritocracy ~ Gini + (1 | YEAR) + (1 | School_Code), data = df)
mod3 <- lm(Belief_in_Meritocracy ~ Ratio_8020, data = df)
mod4 <- lm(Belief_in_Meritocracy ~ pless25zip * pmore75zip, data = df)
export_summs(mod1, mod2, mod3, mod4, digits = 4, error_format = "({std.error})")

#-----------------------------------------------------------------------------
# Table S12.
#-----------------------------------------------------------------------------
df %>%
  filter(YEAR == 2005 | YEAR == 2006 | YEAR == 2007 | YEAR == 2008) %>%
  filter(!is.na(Belief_in_Meritocracy), !is.na(Gini)) %>%
  filter(complete.cases(`Zip Education`, `Zip Income`, `Zip Unemployment`, `Zip White`, `Zip Density`)) %>%
  filter(complete.cases(`Percent Voted Bush`)) %>%
  filter(complete.cases(Female,  White, Citizen, Political_Orientation, `No Religious Affiliation`, Attend_Church,
                        Income, `Mothers Education`, `Fathers Education`, `Conservative Parent Career`)) -> df_limited

q1 <- quantile(df_limited$Gini, 0.01, na.rm = T)
q99 <- quantile(df_limited$Gini, 0.99, na.rm = T)

df_limited %>%
  filter(as.numeric(as.character(Gini)) > q1 & as.numeric(as.character(Gini)) < q99) -> nooutliers

mod1 <- lmer(Belief_in_Meritocracy ~ Gini + `Zip Education` + `Zip Income` + `Zip Unemployment` + `Zip White` + `Zip Density` +
               Female +   White + Citizen + Political_Orientation + `No Religious Affiliation` + Attend_Church +
               Income + `Mothers Education` + `Fathers Education` + `Conservative Parent Career` + `Percent Voted Bush` + 
               (1 | YEAR) + (1 | School_Code), data = nooutliers)

q5 <- quantile(df_limited$Gini, 0.05, na.rm = T)
q95 <- quantile(df_limited$Gini, 0.95, na.rm = T)

df_limited %>%
  filter(as.numeric(as.character(Gini)) > q5 & as.numeric(as.character(Gini)) < q95) -> nooutliers

mod2 <- lmer(Belief_in_Meritocracy ~ Gini + `Zip Education` + `Zip Income` + `Zip Unemployment` + `Zip White` + `Zip Density` +
               Female +   White + Citizen + Political_Orientation + `No Religious Affiliation` + Attend_Church +
               Income + `Mothers Education` + `Fathers Education` + `Conservative Parent Career` + `Percent Voted Bush` + 
               (1 | YEAR) + (1 | School_Code), data = nooutliers)

df %>%
  filter(YEAR == 2005 | YEAR == 2006 | YEAR == 2007 | YEAR == 2008) %>%
  filter(!is.na(Belief_in_Meritocracy), !is.na(Ratio_8020)) %>%
  filter(complete.cases(`Zip Education`, `Zip Income`, `Zip Unemployment`, `Zip White`, `Zip Density`)) %>%
  filter(complete.cases(`Percent Voted Bush`)) %>%
  filter(complete.cases(Female,  White, Citizen, Political_Orientation, `No Religious Affiliation`, Attend_Church,
                        Income, `Mothers Education`, `Fathers Education`, `Conservative Parent Career`)) -> df_limited

q1 <- quantile(df_limited$Ratio_8020, 0.01, na.rm = T)
q99 <- quantile(df_limited$Ratio_8020, 0.99, na.rm = T)

df_limited %>%
  filter(as.numeric(as.character(Ratio_8020)) > q1 & as.numeric(as.character(Ratio_8020)) < q99) -> nooutliers

mod3 <- lmer(Belief_in_Meritocracy ~ Ratio_8020 + `Zip Education` + `Zip Income` + `Zip Unemployment` + `Zip White` + `Zip Density` +
               Female +   White + Citizen + Political_Orientation + `No Religious Affiliation` + Attend_Church +
               Income + `Mothers Education` + `Fathers Education` + `Conservative Parent Career` + `Percent Voted Bush` + 
               (1 | YEAR) + (1 | School_Code), data = nooutliers)

q5 <- quantile(df_limited$Ratio_8020, 0.05, na.rm = T)
q95 <- quantile(df_limited$Ratio_8020, 0.95, na.rm = T)

df_limited %>%
  filter(as.numeric(as.character(Ratio_8020)) > q5 & as.numeric(as.character(Ratio_8020)) < q95) -> nooutliers

mod4 <- lmer(Belief_in_Meritocracy ~ Ratio_8020 + `Zip Education` + `Zip Income` + `Zip Unemployment` + `Zip White` + `Zip Density` +
               Female +   White + Citizen + Political_Orientation + `No Religious Affiliation` + Attend_Church +
               Income + `Mothers Education` + `Fathers Education` + `Conservative Parent Career` + `Percent Voted Bush` + 
               (1 | YEAR) + (1 | School_Code), data = nooutliers)

jtools::export_summs(mod1,mod2,mod3,mod4, digits = 4)

#-----------------------------------------------------------------------------
# Table S13.
#-----------------------------------------------------------------------------
df_clean <- df %>% 
  janitor::clean_names()

mod1 <- lmer(belief_in_meritocracy ~ gini + zip_education + zip_income +
    zip_unemployment + zip_white + zip_density +
    female + white + citizen + political_orientation +
    no_religious_affiliation + attend_church + income +
    mothers_education + fathers_education +
    conservative_parent_career + percent_voted_bush +
    (1 | year) + (1 | school_code),
  data = df_clean, REML = TRUE)

standardize_parameters(mod1, method = "refit") %>%
  as.data.frame() %>%
  mutate(Std_Coefficient = round(Std_Coefficient, digits = 5),
         CI_low = round(CI_low, digits = 4),
         CI_high = round(CI_high, digits = 4)) %>%
  select(-CI)

mod2 <- lmer(belief_in_meritocracy ~ ratio_8020 + zip_education + zip_income +
    zip_unemployment + zip_white + zip_density +
    female + white + citizen + political_orientation +
    no_religious_affiliation + attend_church + income +
    mothers_education + fathers_education +
    conservative_parent_career + percent_voted_bush +
    (1 | year) + (1 | school_code),
  data = df_clean, REML = TRUE
)

standardize_parameters(mod2, method = "refit") %>%
  as.data.frame() %>%
  mutate(Std_Coefficient = round(Std_Coefficient, digits = 5),
         CI_low = round(CI_low, digits = 4),
         CI_high = round(CI_high, digits = 4)) %>%
  select(-CI)

#-----------------------------------------------------------------------------
# Figure S1.
#-----------------------------------------------------------------------------
interaction_data_std <- df %>%
  mutate(Zip_Education = `Zip Education`,
    Zip_Income = `Zip Income`,
    Zip_Unemployment = `Zip Unemployment`,
    Zip_White = `Zip White`,
    Zip_Density = `Zip Density`,
    Percent_Voted_Bush = `Percent Voted Bush`,
    School_Code = `School Code`,
    No_Religion = `No Religious Affiliation`,
    Mothers_Education = `Mothers Education`,
    Fathers_Education = `Fathers Education`,
    Conservative_Parent = `Conservative Parent Career`) %>%
  select(YEAR, Income, School_Code, Belief_in_Meritocracy, Gini, 
         Zip_Education, Zip_Income, Zip_Unemployment, Zip_White, Zip_Density,
         Percent_Voted_Bush, No_Religion, 
         Mothers_Education, Fathers_Education, Conservative_Parent,
         Female, Age, White, Citizen, Political_Orientation, Attend_Church, HOMEZIP, School_Code, YEAR,
         pmore75zip, pless25zip) %>%
  filter(complete.cases(Belief_in_Meritocracy, Gini, Income)) %>%
  mutate(across(c(Belief_in_Meritocracy, Gini, Income,
      Zip_Education, Zip_Income, Zip_Unemployment, Zip_White, Zip_Density,
      Percent_Voted_Bush, Political_Orientation, Attend_Church, No_Religion,
      Mothers_Education, Fathers_Education, Conservative_Parent, Age, pless25zip),
    ~ scale(.)[, 1])) %>%
  as.data.frame()

PlotS1 <- interflex(
  estimator = "binning",
  Y = "Belief_in_Meritocracy",
  D = "pless25zip",
  X = "pmore75zip",
  Z = c("Female", "White", "Citizen",
        "Political_Orientation", "Attend_Church", "No_Religion",
        "Mothers_Education", "Fathers_Education", "Conservative_Parent",
        "Zip_Education", "Zip_Income", "Zip_Unemployment", "Zip_White",
        "Zip_Density", "Percent_Voted_Bush"),
  cl = c("YEAR", "School_Code"),
  data = interaction_data_std, na.rm = TRUE,
  Ylabel = "", Dlabel = "GINI", xlab = "",
  cex.main = 1.2,
  theme.bw = TRUE, show.grid = FALSE,
  parallel=TRUE)

data <- PlotS1$est.bin[["pless25zip=-0.22 (50%)"]] %>%
  as_data_frame()

ggplot(data, aes(x = x0, y = coef)) +
  geom_point(color = "red", size = 3, fill = "white") + 
  geom_errorbar(aes(ymin = CI.lower, ymax = CI.upper), width = 0.1, color = "red", size = 0.75) +
  geom_hline(yintercept = 0, linetype = "solid", color = "grey") +
  scale_x_continuous(limits = c(0, 1), breaks = seq(0, 1, by = 0.25)) +
  scale_y_continuous(limits = c(-0.2, 0.01), breaks = seq(-0.2, 0.01, by = 0.1)) + 
  theme_bw() + 
  labs(x = "", y = "Marginal effect of % Below $25k") +
  theme(axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 15),
        axis.title.x = element_text(size = 15),
        axis.title.y = element_text(size = 15),
        panel.border = element_blank()) +
  annotation_custom(ggplotGrob(
      ggplot(interaction_data_std, aes(x = Income)) +
        geom_histogram(bins = 100, fill = "grey", color = "black") +
        theme_void() + 
        theme(plot.margin = margin(t = -1, b = -1))),
    xmin = 0, xmax = 1, ymin = -0.2, ymax = -0.15)

#-----------------------------------------------------------------------------
# Table S14. Standardized coefficient estimates from the Interflex binning model above in Figure S1.
#-----------------------------------------------------------------------------
round(summary(PlotS1$model.binnin)$coefficients, digits = 5)

#-----------------------------------------------------------------------------
# Figure S2.
#-----------------------------------------------------------------------------
interaction_data_std <- df %>%
  mutate(Zip_Education = `Zip Education`,
    Zip_Income = `Zip Income`,
    Zip_Unemployment = `Zip Unemployment`,
    Zip_White = `Zip White`,
    Zip_Density = `Zip Density`,
    Percent_Voted_Bush = `Percent Voted Bush`,
    School_Code = `School Code`,
    No_Religion = `No Religious Affiliation`,
    Mothers_Education = `Mothers Education`,
    Fathers_Education = `Fathers Education`,
    Conservative_Parent = `Conservative Parent Career`) %>%
  select(YEAR, Income, School_Code, Belief_in_Meritocracy, Gini, 
         Zip_Education, Zip_Income, Zip_Unemployment, Zip_White, Zip_Density,
         Percent_Voted_Bush, No_Religion, 
         Mothers_Education, Fathers_Education, Conservative_Parent,
         Female, Age, White, Citizen, Political_Orientation, Attend_Church, HOMEZIP, School_Code, YEAR) %>%
  filter(complete.cases(Belief_in_Meritocracy, Gini, Income)) %>%
  mutate(across(c(Belief_in_Meritocracy, Gini,
      Zip_Education, Zip_Income, Zip_Unemployment, Zip_White, Zip_Density,
      Percent_Voted_Bush, Political_Orientation, Attend_Church, No_Religion,
      Mothers_Education, Fathers_Education, Conservative_Parent, Age),
    ~ scale(.)[, 1])) %>%
  as.data.frame()

PlotS2 <- interflex(
  estimator = "binning",
  Y = "Belief_in_Meritocracy",
  D = "Gini",
  X = "Income",
  Z = c("Female", "White", "Citizen",
        "Political_Orientation", "Attend_Church", "No_Religion",
        "Mothers_Education", "Fathers_Education", "Conservative_Parent",
        "Zip_Education", "Zip_Income", "Zip_Unemployment", "Zip_White",
        "Zip_Density", "Percent_Voted_Bush"),
  cl = c("YEAR", "School_Code"),
  data = interaction_data_std, na.rm = TRUE,
  Ylabel = "Belief in Meritocracy", Dlabel = "GINI", xlab = "",
  cex.main = 1.2,
  theme.bw = TRUE, show.grid = FALSE,
  parallel=TRUE)

data <- PlotS2$est.bin[["Gini=-0.083 (50%)"]] %>%
  as_data_frame()

ggplot(data, aes(x = x0, y = coef)) +
  geom_point(color = "red", size = 3, fill = "white") +
  geom_errorbar(aes(ymin = CI.lower, ymax = CI.upper), width = 0.1, color = "red", size = 0.75) + 
  geom_hline(yintercept = 0, linetype = "solid", color = "grey") + 
  scale_x_continuous(limits = c(0, 1), breaks = seq(0, 1, by = 0.25)) +
  scale_y_continuous(limits = c(-0.05, 0.01), breaks = seq(-0.05, 0.01, by = 0.1)) + 
  theme_bw() + 
  labs(x = "", y = "Marginal effect of Gini Coefficient") +
  theme(axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 15),
        axis.title.x = element_text(size = 15),
        axis.title.y = element_text(size = 15),
        panel.border = element_blank()) +
  annotation_custom(ggplotGrob(ggplot(interaction_data_std, aes(x = Income)) +
        geom_histogram(bins = 100, fill = "grey", color = "black") + 
        theme_void() + 
        theme(plot.margin = margin(t = -1, b = -1))),
    xmin = 0, xmax = 1, ymin = -0.05, ymax = -0.043)

#-----------------------------------------------------------------------------
# Table S15. Standardized coefficient estimates from the Interflex binning model above in Figure S2.
#-----------------------------------------------------------------------------
round(summary(PlotS2$model.binning)$coefficients, digits = 5)

#-----------------------------------------------------------------------------
# Figure S3.
#-----------------------------------------------------------------------------
interaction_data_std <- df %>%
  mutate(Zip_Education = `Zip Education`,
    Zip_Income = `Zip Income`,
    Zip_Unemployment = `Zip Unemployment`,
    Zip_White = `Zip White`,
    Zip_Density = `Zip Density`,
    Percent_Voted_Bush = `Percent Voted Bush`,
    School_Code = `School Code`,
    No_Religion = `No Religious Affiliation`,
    Mothers_Education = `Mothers Education`,
    Fathers_Education = `Fathers Education`,
    Conservative_Parent = `Conservative Parent Career`) %>%
  select(YEAR, Income, School_Code, Belief_in_Meritocracy, Ratio_8020, 
         Zip_Education, Zip_Income, Zip_Unemployment, Zip_White, Zip_Density,
         Percent_Voted_Bush, No_Religion, 
         Mothers_Education, Fathers_Education, Conservative_Parent,
         Female, Age, White, Citizen, Political_Orientation, Attend_Church, HOMEZIP, School_Code, YEAR) %>%
  filter(complete.cases(Belief_in_Meritocracy, Ratio_8020, Income)) %>%
  mutate(across(c(Belief_in_Meritocracy, Ratio_8020,
      Zip_Education, Zip_Income, Zip_Unemployment, Zip_White, Zip_Density,
      Percent_Voted_Bush, Political_Orientation, Attend_Church, No_Religion,
      Mothers_Education, Fathers_Education, Conservative_Parent, Age),
    ~ scale(.)[, 1])) %>%
  as.data.frame()

PlotS3 <- interflex(
  estimator = "binning",
  Y = "Belief_in_Meritocracy",
  D = "Ratio_8020",
  X = "Income",
  Z = c("Female", "White", "Citizen",
        "Political_Orientation", "Attend_Church", "No_Religion",
        "Mothers_Education", "Fathers_Education", "Conservative_Parent",
        "Zip_Education", "Zip_Income", "Zip_Unemployment", "Zip_White",
        "Zip_Density", "Percent_Voted_Bush"),
  cl = c("YEAR", "School_Code"),
  data = interaction_data_std, na.rm = TRUE,
  Ylabel = "Belief in Meritocracy", Dlabel = "80:20", xlab = "",
  cex.main = 1.2,
  theme.bw = TRUE, show.grid = FALSE,
  parallel=TRUE)

data <- PlotS3$est.bin[["Ratio_8020=-0.087 (50%)"]] %>%
  as_data_frame()

ggplot(data, aes(x = x0, y = coef)) +
  geom_point(color = "red", size = 3, fill = "white") + 
  geom_errorbar(aes(ymin = CI.lower, ymax = CI.upper), width = 0.1, color = "red", size = 0.75) + 
  geom_hline(yintercept = 0, linetype = "solid", color = "grey") +
  scale_x_continuous(limits = c(0, 1), breaks = seq(0, 1, by = 0.25)) +
  scale_y_continuous(limits = c(-0.05, 0.01), breaks = seq(-0.05, 0.01, by = 0.1)) + 
  theme_bw() + 
  labs(x = "", y = "Marginal effect of 80:20") +
  theme(axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 15),
        axis.title.x = element_text(size = 15),
        axis.title.y = element_text(size = 15),
        panel.border = element_blank()) +
  annotation_custom(ggplotGrob(ggplot(interaction_data_std, aes(x = Income)) +
        geom_histogram(bins = 100, fill = "grey", color = "black") +
        theme_void() + 
        theme(plot.margin = margin(t = -1, b = -1))),
    xmin = 0, xmax = 1, ymin = -0.05, ymax = -0.043)

#-----------------------------------------------------------------------------
# Table S16. Standardized coefficient estimates from the Interflex binning model above in Figure S3.
#-----------------------------------------------------------------------------
round(summary(PlotS3$model.binning)$coefficients, digits = 5)

#-----------------------------------------------------------------------------
# Table S17.
#-----------------------------------------------------------------------------
mod1 <- lmer(Belief_in_Meritocracy ~ Gini + `Zip Education` + `Zip Income` + `Zip Unemployment` + `Zip White` + `Zip Density` +
                Female +   White + Citizen + Political_Orientation + `No Religious Affiliation` + Attend_Church +
                Income + `Mothers Education` + `Fathers Education` + `Conservative Parent Career` + `Percent Voted Bush` + Selectivity +
                (1 | YEAR) + (1 | School_Code), data = df)
mod2 <-  lmer(Belief_in_Meritocracy ~ Ratio_8020 + `Zip Education` + `Zip Income` + `Zip Unemployment` + `Zip White` + `Zip Density` +
                 Female +   White + Citizen + Political_Orientation + `No Religious Affiliation` + Attend_Church +
                 Income + `Mothers Education` + `Fathers Education` + `Conservative Parent Career` + `Percent Voted Bush` + Selectivity +
                 (1 | YEAR) + (1 | School_Code), data = df)
mod3 <-  lmer(Belief_in_Meritocracy ~ pless25zip * pmore75zip + `Zip Education` + `Zip Income` + `Zip Unemployment` + `Zip White` + `Zip Density` +
                 Female +   White + Citizen + Political_Orientation + `No Religious Affiliation` + Attend_Church +
                 Income + `Mothers Education` + `Fathers Education` + `Conservative Parent Career` + `Percent Voted Bush` + Selectivity +
                 (1 | YEAR) + (1 | School_Code), data = df)

export_summs(mod1, mod2, mod3, digits = 4, error_format = "({std.error})")

# ---------------------------------------------------------------------------------------------------------------
# ------------------------ Appendix-E: Alternative Interaction Model Results  -----------------------------------
# ---------------------------------------------------------------------------------------------------------------

#-----------------------------------------------------------------------------
# Table S18.
#-----------------------------------------------------------------------------
df %>%
  mutate(Zip_Education = `Zip Education`,
    Zip_Income = `Zip Income`,
    Zip_Unemployment = `Zip Unemployment`,
    Zip_White = `Zip White`,
    Zip_Density = `Zip Density`,
    Percent_Voted_Bush = `Percent Voted Bush`,
    School_Code = `School Code`,
    No_Religion = `No Religious Affiliation`,
    Mothers_Education = `Mothers Education`,
    Fathers_Education = `Fathers Education`,
    Conservative_Parent = `Conservative Parent Career`) %>%
  select(YEAR, Income, School_Code, Belief_in_Meritocracy, Gini, 
         Zip_Education, Zip_Income, Zip_Unemployment, Zip_White, Zip_Density,
         Percent_Voted_Bush, No_Religion, 
         Mothers_Education, Fathers_Education, Conservative_Parent,
         Female, Age, White, Citizen, Political_Orientation, Attend_Church, Ratio_8020) %>%
  filter(complete.cases(Belief_in_Meritocracy, Ratio_8020, Income)) %>%
  as.data.frame() -> interaction_data

interflex(estimator = "binning",Y = "Belief_in_Meritocracy", D = "Gini", X = "Income", 
          Z = c("Female", "White", "Citizen", "Income",
                "Political_Orientation", "Attend_Church", "No_Religion",
                "Mothers_Education", "Fathers_Education", "Conservative_Parent",
                "Zip_Education", "Zip_Income", "Zip_Unemployment", "Zip_White",
                "Zip_Density", "Percent_Voted_Bush"), cl = c("YEAR", "School_Code"),
          data = interaction_data, na.rm = TRUE,
          Ylabel = "Believe Meritocracy", Dlabel = "Gini", Xlabel="Parent's Income",
          main = "Marginal Effects", cex.main = 1.2, nbins = 4,
          theme.bw = TRUE, show.grid = FALSE,
          parallel=TRUE)$est.bin

#-----------------------------------------------------------------------------
# Table S19.
#-----------------------------------------------------------------------------
interflex(estimator = "binning",Y = "Belief_in_Meritocracy", D = "Gini", X = "Income", 
          Z = c("Female", "White", "Citizen", "Income",
                "Political_Orientation", "Attend_Church", "No_Religion",
                "Mothers_Education", "Fathers_Education", "Conservative_Parent",
                "Zip_Education", "Zip_Income", "Zip_Unemployment", "Zip_White",
                "Zip_Density", "Percent_Voted_Bush"), cl = c("YEAR", "School_Code"),
          data = interaction_data, na.rm = TRUE,
          Ylabel = "Believe Meritocracy", Dlabel = "Gini", Xlabel="Parent's Income",
          main = "Marginal Effects", cex.main = 1.2, nbins = 5,
          theme.bw = TRUE, show.grid = FALSE,
          parallel=TRUE)$est.bin

#-----------------------------------------------------------------------------
# Table S20.
#-----------------------------------------------------------------------------
interflex(estimator = "binning",Y = "Belief_in_Meritocracy", D = "Ratio_8020", X = "Income", 
                   Z = c("Female", "White", "Citizen", "Income",
                         "Political_Orientation", "Attend_Church", "No_Religion",
                         "Mothers_Education", "Fathers_Education", "Conservative_Parent",
                         "Zip_Education", "Zip_Income", "Zip_Unemployment", "Zip_White",
                         "Zip_Density", "Percent_Voted_Bush"), cl = c("YEAR", "School_Code"),
                   data = interaction_data, na.rm = TRUE,
                   Ylabel = "Believe Meritocracy", Dlabel = "80/20 Ratio", Xlabel="Parent's Income",
                   main = "Marginal Effects", cex.main = 1.2, 
                   theme.bw = TRUE, show.grid = FALSE,
                   parallel=TRUE)$est.bin

#-----------------------------------------------------------------------------
# Table S21.
#-----------------------------------------------------------------------------
interflex(estimator = "binning",Y = "Belief_in_Meritocracy", D = "Ratio_8020", X = "Income", 
                   Z = c("Female", "White", "Citizen", "Income",
                         "Political_Orientation", "Attend_Church", "No_Religion",
                         "Mothers_Education", "Fathers_Education", "Conservative_Parent",
                         "Zip_Education", "Zip_Income", "Zip_Unemployment", "Zip_White",
                         "Zip_Density", "Percent_Voted_Bush"), cl = c("YEAR", "School_Code"),
                   data = interaction_data, na.rm = TRUE,
                   Ylabel = "Believe Meritocracy", Dlabel = "80/20 Ratio", Xlabel="Parent's Income",
                   main = "Marginal Effects", cex.main = 1.2, nbins = 4,
                   theme.bw = TRUE, show.grid = FALSE,
                   parallel=TRUE)$est.bin
#-----------------------------------------------------------------------------
# Table S22.
#-----------------------------------------------------------------------------
interflex(estimator = "binning",Y = "Belief_in_Meritocracy", D = "Ratio_8020", X = "Income", 
                   Z = c("Female", "White", "Citizen", "Income",
                         "Political_Orientation", "Attend_Church", "No_Religion",
                         "Mothers_Education", "Fathers_Education", "Conservative_Parent",
                         "Zip_Education", "Zip_Income", "Zip_Unemployment", "Zip_White",
                         "Zip_Density", "Percent_Voted_Bush"), cl = c("YEAR", "School_Code"),
                   data = interaction_data, na.rm = TRUE,
                   Ylabel = "Believe Meritocracy", Dlabel = "80/20 Ratio", Xlabel="Parent's Income",
                   main = "Marginal Effects", cex.main = 1.2, nbins = 5,
                   theme.bw = TRUE, show.grid = FALSE,
                   parallel=TRUE)$est.bin

# ---------------------------------------------------------------------------------------------------------------
# ------------------------   Appendix-G: Empathy as a Potential Moderator.    -----------------------------------
# ---------------------------------------------------------------------------------------------------------------

#-----------------------------------------------------------------------------
# Figure S4.
#-----------------------------------------------------------------------------
df %>%
  filter(UnderstandingofOthers == 0 | UnderstandingofOthers == 0.25) -> LowEmpathy

df %>%
  filter(UnderstandingofOthers == 1 | UnderstandingofOthers == 0.75) -> HighEmpathy

LowEmpathy %>%
  mutate(Gini = Gini,
    Zip_Education = `Zip Education`,
    Zip_Income = `Zip Income`,
    Zip_Unemployment = `Zip Unemployment`,
    Zip_White = `Zip White`,
    Zip_Density = `Zip Density`,
    Percent_Voted_Bush = `Percent Voted Bush`,
    School_Code = `School Code`,
    No_Religion = `No Religious Affiliation`,
    Mothers_Education = `Mothers Education`,
    Fathers_Education = `Fathers Education`,
    Conservative_Parent = `Conservative Parent Career`) %>%
  select(YEAR, Income, School_Code, Belief_in_Meritocracy, Gini, Ratio_8020,
         Zip_Education, Zip_Income, Zip_Unemployment, Zip_White, Zip_Density,
         Percent_Voted_Bush, No_Religion, 
         Mothers_Education, Fathers_Education, Conservative_Parent,
         Female, Age, White, Citizen, Political_Orientation, Attend_Church) %>%
  filter(complete.cases(Belief_in_Meritocracy, Gini, Income)) %>%
  as.data.frame() -> interaction_data1

PlotS4A <- interflex(estimator = "binning",Y = "Belief_in_Meritocracy", D = "Gini", X = "Income", 
                   Z = c("Female", "White", "Citizen",
                         "Political_Orientation", "Attend_Church", "No_Religion",
                         "Mothers_Education", "Fathers_Education", "Conservative_Parent",
                         "Zip_Education", "Zip_Income", "Zip_Unemployment", "Zip_White",
                         "Zip_Density", "Percent_Voted_Bush"),cl = c("YEAR", "School_Code"),
                   data = interaction_data1, na.rm = TRUE,
                   ylab = "", Dlabel = "GINI", xlab = "", ylim = c(-0.4, 0.2),
                   cex.main = 1.2, 
                   theme.bw = TRUE, show.grid = FALSE,
                   parallel=TRUE)

data <- PlotS4A$est.bin[["Gini=0.594 (50%)"]] %>%
  as_data_frame()

ggplot(data, aes(x = x0, y = coef)) +
  geom_point(color = "red", size = 3, fill = "white") +
  geom_errorbar(aes(ymin = CI.lower, ymax = CI.upper), width = 0.1, color = "red", size = 0.75) +
  geom_hline(yintercept = 0, linetype = "solid", color = "grey") +
  scale_x_continuous(limits = c(0, 1), breaks = seq(0, 1, by = 0.25)) +
  scale_y_continuous(limits = c(-0.8, 0.2), breaks = seq(-0.8, 0.2, by = 0.1)) + 
  theme_bw() + 
  labs(x = "Low Empathy", y = "Moderator: Parental Income") +
  theme(axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 15),
        axis.title.x = element_text(size = 15),
        axis.title.y = element_text(size = 15),
        panel.border = element_blank()) +
  annotation_custom(
    ggplotGrob(
      ggplot(interaction_data1, aes(x = Income)) +
        geom_histogram(bins = 100, fill = "grey", color = "black") +
        theme_void() +
        theme(plot.margin = margin(t = -1, b = -1))
    ),
    xmin = 0, xmax = 1, ymin = -0.8, ymax = -0.6 
  ) -> PlotS4A

HighEmpathy %>%
  mutate(Gini = Gini,
    Zip_Education = `Zip Education`,
    Zip_Income = `Zip Income`,
    Zip_Unemployment = `Zip Unemployment`,
    Zip_White = `Zip White`,
    Zip_Density = `Zip Density`,
    Percent_Voted_Bush = `Percent Voted Bush`,
    School_Code = `School Code`,
    No_Religion = `No Religious Affiliation`,
    Mothers_Education = `Mothers Education`,
    Fathers_Education = `Fathers Education`,
    Conservative_Parent = `Conservative Parent Career`) %>%
  select(YEAR, Income, School_Code, Belief_in_Meritocracy, Gini, Ratio_8020,
         Zip_Education, Zip_Income, Zip_Unemployment, Zip_White, Zip_Density,
         Percent_Voted_Bush, No_Religion, 
         Mothers_Education, Fathers_Education, Conservative_Parent,
         Female, Age, White, Citizen, Political_Orientation, Attend_Church) %>%
  filter(complete.cases(Belief_in_Meritocracy, Gini, Income)) %>%
  as.data.frame() -> interaction_data2

PlotS4B <- interflex(estimator = "binning",Y = "Belief_in_Meritocracy", D = "Gini", X = "Income", 
                   Z = c("Female", "White", "Citizen",
                         "Political_Orientation", "Attend_Church", "No_Religion",
                         "Mothers_Education", "Fathers_Education", "Conservative_Parent",
                         "Zip_Education", "Zip_Income", "Zip_Unemployment", "Zip_White",
                         "Zip_Density", "Percent_Voted_Bush"),cl = c("YEAR", "School_Code"),
                   data = interaction_data2, na.rm = TRUE,
                   ylab = "", Dlabel = "GINI", xlab = "", ylim = c(-0.4, 0.2),
                   cex.main = 1.2, 
                   theme.bw = TRUE, show.grid = FALSE,
                   parallel=TRUE)

data <- PlotS4B$est.bin[["Gini=0.591 (50%)"]] %>%
  as_data_frame()

ggplot(data, aes(x = x0, y = coef)) +
  geom_point(color = "red", size = 3, fill = "white") +
  geom_errorbar(aes(ymin = CI.lower, ymax = CI.upper), width = 0.1, color = "red", size = 0.75) +
  geom_hline(yintercept = 0, linetype = "solid", color = "grey") +
  scale_x_continuous(limits = c(0, 1), breaks = seq(0, 1, by = 0.25)) +
  scale_y_continuous(limits = c(-0.8, 0.2), breaks = seq(-0.8, 0.2, by = 0.1)) + 
  theme_bw() + 
  labs(x = "Low Empathy", y = "Moderator: Parental Income") +
  theme(axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 15),
        axis.title.x = element_text(size = 15),
        axis.title.y = element_text(size = 15),
        panel.border = element_blank()) +
  annotation_custom(
    ggplotGrob(
      ggplot(interaction_data2, aes(x = Income)) +
        geom_histogram(bins = 100, fill = "grey", color = "black") +
        theme_void() +
        theme(plot.margin = margin(t = -1, b = -1))
    ),
    xmin = 0, xmax = 1, ymin = -0.8, ymax = -0.6 
  ) -> PlotS4B


ggarrange(PlotS4A, PlotS4B)

#-----------------------------------------------------------------------------
# Figure S5.
#-----------------------------------------------------------------------------
PlotS5A <- interflex(estimator = "binning",Y = "Belief_in_Meritocracy", D = "Ratio_8020", X = "Income", 
                   Z = c("Female", "White", "Citizen",
                         "Political_Orientation", "Attend_Church", "No_Religion",
                         "Mothers_Education", "Fathers_Education", "Conservative_Parent",
                         "Zip_Education", "Zip_Income", "Zip_Unemployment", "Zip_White",
                         "Zip_Density", "Percent_Voted_Bush"),cl = c("YEAR", "School_Code"),
                   data = interaction_data1, na.rm = TRUE,
                   ylab = "", Dlabel = "GINI", xlab = "", ylim = c(-0.4, 0.2),
                   cex.main = 1.2, 
                   theme.bw = TRUE, show.grid = FALSE,
                   parallel=TRUE)

data <- PlotS5A$est.bin[["Ratio_8020=0.073 (50%)"]] %>%
  as_data_frame()

ggplot(data, aes(x = x0, y = coef)) +
  geom_point(color = "red", size = 3, fill = "white") +
  geom_errorbar(aes(ymin = CI.lower, ymax = CI.upper), width = 0.1, color = "red", size = 0.75) +
  geom_hline(yintercept = 0, linetype = "solid", color = "grey") +
  scale_x_continuous(limits = c(0, 1), breaks = seq(0, 1, by = 0.25)) +
  scale_y_continuous(limits = c(-0.8, 0.2), breaks = seq(-0.8, 0.2, by = 0.1)) + 
  theme_bw() + 
  labs(x = "Low Empathy", y = "Marginal effect of Parental Income") +
  theme(axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 15),
        axis.title.x = element_text(size = 15),
        axis.title.y = element_text(size = 15),
        panel.border = element_blank()) +
  annotation_custom(ggplotGrob(
      ggplot(interaction_data1, aes(x = Income)) +
        geom_histogram(bins = 100, fill = "grey", color = "black") +
        theme_void() +
        theme(plot.margin = margin(t = -1, b = -1))),
    xmin = 0, xmax = 1, ymin = -0.8, ymax = -0.6) -> PlotS5A

PlotS5B <- interflex(estimator = "binning",Y = "Belief_in_Meritocracy", D = "Ratio_8020", X = "Income", 
                   Z = c("Female", "White", "Citizen",
                         "Political_Orientation", "Attend_Church", "No_Religion",
                         "Mothers_Education", "Fathers_Education", "Conservative_Parent",
                         "Zip_Education", "Zip_Income", "Zip_Unemployment", "Zip_White",
                         "Zip_Density", "Percent_Voted_Bush"),cl = c("YEAR", "School_Code"),
                   data = interaction_data2, na.rm = TRUE,
                   ylab = "", Dlabel = "GINI", xlab = "", ylim = c(-0.4, 0.2),
                   cex.main = 1.2, 
                   theme.bw = TRUE, show.grid = FALSE,
                   parallel=TRUE)

data <- PlotS5B$est.bin[["Ratio_8020=0.073 (50%)"]] %>%
  as_data_frame()

ggplot(data, aes(x = x0, y = coef)) +
  geom_point(color = "red", size = 3, fill = "white") +
  geom_errorbar(aes(ymin = CI.lower, ymax = CI.upper), width = 0.1, color = "red", size = 0.75) +
  geom_hline(yintercept = 0, linetype = "solid", color = "grey") +
  scale_x_continuous(limits = c(0, 1), breaks = seq(0, 1, by = 0.25)) +
  scale_y_continuous(limits = c(-0.8, 0.2), breaks = seq(-0.8, 0.2, by = 0.1)) + 
  theme_bw() + 
  labs(x = "High Empathy", y = "") +
  theme(axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 15),
        axis.title.x = element_text(size = 15),
        axis.title.y = element_text(size = 15),
        panel.border = element_blank()) +
  annotation_custom(ggplotGrob(
      ggplot(interaction_data2, aes(x = Income)) +
        geom_histogram(bins = 100, fill = "grey", color = "black") +
        theme_void() +
        theme(plot.margin = margin(t = -1, b = -1))
    ),
    xmin = 0, xmax = 1, ymin = -0.8, ymax = -0.6
  ) -> PlotS5B; PlotS5B

ggarrange(PlotS5A, PlotS5B)
